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Decomposition algorithms 
for globally solving mathematical programs 

with affine equilibrium constraints u 

L. D. MuiS • T. D. Quocl • L. T. H. Aifl • P. D. Tao| 

o . 

Abstract A mathematical programming problem with affine equilibrium constraints (AM- 
PEC) is a bilevel programming problem where the lower one is a parametric affine varia- 
tional inequality. We formulate some classes of bilevel programming in forms of AMPEC. 
Then we use a regularization technique to formulate the resulting problem as a mathemati- 
cal program with an additional constraint defined by the difference of two convex functions 
(DC function). A main feature of this DC decomposition is that the second component de- 
pends upon only the parameter in the lower problem. This property allows us to develop 
branch-and-bound algorithms for globally solving AMPEC where the adaptive rectangular 
bisection takes place only in the space of the parameter. As an example, we use the pro- 
posed algorithm to solve a bilevel Nash-Cournot equilibrium market model. Computational 
results show the efficiency of the proposed algorithm. 

Keywords Mathematical programs with affine equilibrium constraints ■ regularization 
• bilevel convex quadratic programming ■ DC formulation ■ global optimization ■ Nash- 
^ ! Cournot model. 

m 
cn 

]Q 1 Introduction 

o ■ 

We consider the following mathematical programming problem with affine (not necessarily 
monotone) variational inequality constraints, that we call shortly AMPEC: 

min f(x,y) (P) 

s.t. (x,y)eS, (1) 
x E C, (Ax + By + a) T (v - x) > 0, \fv G C, (2) 

where R n+m , ^ C C W n are two closed convex sets, / : R m+n -> R 

is a convex function, A, B are given appropriate real matrices, and a E R n . This class 
of optimization problems is known to be very difficult to solve due to its nonconvexity, 
nondifferentiability and loss of constraint qualification. However such problems arise fre- 
quently in applications, for example, in shape optimization, design transportation network, 
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economic modeling and data mining. A natural way to handle a nested problem such as 
Problem © is to reduce it into an one-level optimization problem by using the Krush- 
Kuhn-Tucker theorem for the lower variational inequality. Several algorithms for globally 
solving the reduced mathematical programs with complementarity constraints are proposed 
(see, e.g. (3l[5l[T9l|22l). Since the number of the complementarity constraints is just equal 
to the number of constraints defining the set C in the lower variational inequality problem, 
these global optimization algorithms become expensive when the number of constraints is 
high, for example, when C := {x E R n | x > 0, Cj(x) < 0, j — 1, . . . ,p} (often appears 
in practice) with either n or p are somewhat large. 

In this paper, we propose another solution-approach to AMPEC without using the 
Krush-Kuhn-Tucker theorem for the lower variational inequality. Instead, we use a reg- 
ularization technique to formulate AMEC as a mathematical program with an additional 
constraint defined by gi(x,y) — hi(x,y) < 0, where gi and hi are differentiable con- 
vex functions. The main feature of this constraint is that the second component h\ can 
be chosen such a way so that it depends upon only the parameter y. Moreover, in some 
special important cases such as bilevel convex quadratic problems, hi is separable. This 
formulation allows us to develop decomposition branch-and-bound algorithms for glob- 
ally solving AMPEC where the branching operation involving only the parameter in the 
lower variational inequality. Unlike the existing global optimization algorithms mentioned 
above, the proposed algorithms can solve AMPEC where the constraint set C is given as 
C := {x E R n | x > 0, Cj(x) < 0, j = 1, . . . ,p} with n and p relatively large. As 
an example, we use the proposed algorithm to find a global optimal equilibrium pair to 
a bilevel Nash-Coumot equilibrium market model. We tested the proposed algorithm by 
some randomly generated data. The numerical results show that our algorithm can solve 
this bilevel model with high dimension. 

The paper is organized as follows. In the next section we give DC formulations to AM- 
PEC by using suitable regularization matrices. Some important special cases of AMPEC 
are presented at the end of this section. The third section is devoted to description of a 
branch-and-bound algorithm for globally solving a bilevel Nash-Coumot equilibrium mar- 
ket model by using a DC decomposition, where the second component is separable and 
depends upon only the parameter y. We close the paper with some computational experi- 
ences and results. 

2 DC Formulations and Examples 

In Problem ©, as usual, we will refer to x as a primary variable or decision variable and 
y as a parameter. We call (x, y) a feasible point to © if (x, y) E S and x solves the 
lower variational inequality ©. Note that when A is symmetric positive semidefinite, the 
variational inequality © is equivalent to the parametric convex quadratic problem 

min {(p(x, y) := -x T Ax + (By + a) T x \ x E C}. (3) 
In this case, Problem (10) becomes a bilevel convex program 

mm{f(x,y)\ (x,y) E S}, (BP) 
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where x solves the convex quadratic program 



mm 



{<p{x,y) 



:= -x T Ax + (By + afx \ x G C) . 



(4) 



In the general case, when A is indefinite, the variational inequality (O is not necessarily 
equivalent to the problem ©. So Problem ©, in general, cannot be reformulated as a 
bilevel problem of the form (ED IT71II51. 

2.1 DC Formulations 

The main difficulty of Problem © is that the constraint defined by the variational inequality 
© is neither convex nor given explicitly as a constrained set of an ordinary mathematical 
programming problem. A natural way is to reduce Problem ()P) to an ordinary mathematical 
programming problem. We will reformulate Problem (0 as a smoothly DC program. For 
this purpose, we use a gap function introduced in fl2T1 to formulate the variational inequality 
© as an equation defined by a smoothly DC function. We recall that a function / is said 
to be DC on a convex set D if it can be expressed as the difference of two convex functions 
on D, i.e. f = g — h, where g and h are convex on D. 

More precisely, for each (x,y), following the idea from ETI we define the function 
g(x, y) by setting 



where G is an arbitrary symmetric positive definite (n x n) -matrix. We refer to G as a 
regularization matrix. Since G is positive definite, the problem defining g(x, y) is uniquely 
solvable for every (x, y), i.e. g is well-defined. 

The following lemma gives the properties of the gap function g whose proof can be 
done similarly to the proof of Lemma 2.1 in [1211 . 

Lemma 2.1. Let g be given by ©. Then 

(i) g(x, y) > Ofor every (x, y) G C x R m , 

(ii) (x, y) G S, x G C, g(x, y) = if and only if (x, y) is feasible solution to Problem 



The following proposition shows that, with a suitable choice of the regularization ma- 
trix G, the function g can be decomposed as the difference of two convex functions (DC 
function). Note that any symmetric matrix A can be expressed as A = A\ — A 2 , where Ai 
is symmetric positive definite and A 2 is symmetric. In what follows by diag(ct) we denote 
the diagonal matrix whose every diagonal entry is a. 

Proposition 2.1. Suppose that A is symmetric and A = Ax— A 2 , where Ai is a symmetric 
positive definite matrix and A 2 is a symmetric matrix such that A 2 + \U T U is positive 
(semi)definite, and U, V are two appropriate matrices satisfying U T V = B. Let G = 2A\. 



9{x,y) : 



= max 




(5) 



©. 



Then 



g(x,y) = gi(x,y) - h(x,y), 
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where g± and hi are two differentiable convex functions given by 

9i(x,y) = ^\\Ux + Vy\\ 2 + a T x 

+ max { [(A 1 + A 2 )x -By- afv - v T A 1 v], (7) 

and 

h(x, y) = l -x T {2A 2 + U T U)x + l -\\Vy\\ 2 . (8) 
Proof. With a simple arrangement from ©, it shows that 

g(x, y) = x T Ax — -x T Gx + x 7 By + a T x 

+ max{— v T Ax — v 7 By — a T v v T Gv + x T Gv}. (9) 

vec 2 

Since A = Ai — A 2 and G = 2Ai, the last expression implies 

g{x, y) = —x T A 2 x + x 7 By + a T x 

+ max{-t; T Ai^ + [(A 1 + A 2 )x - By - a] 7 v}. (10) 

v<=C 

On the other hand, since B = U T V we can express 

2x T By = 2x T U T Vy = \\Ux + Vy\\ 2 - \\Ux\\ 2 - \\Vy\\ 2 . 
Substituting this expression into (flOl) we get 

g(x,y) = i\\Ux + Vy\\ 2 - ~\\Ux\\ 2 - l -\\Vy\\ 2 + a T x 
+ max{-w T Ait; + [(A 1 + A 2 )x - By - a} T v}. 

Hence, 

g(x,y) = gi(x,y) - hi(x,y), 

where g x and h\ are two functions given by © and ([8]), respectively. Since A 2 + \U 7 U 
is positive semidefinite, h\ is convex. Clearly, h\ is differentiable everywhere, while g\ is 
differentiable everywhere because the convex program (strongly quadratic concave maxi- 
mization): 

max { - v T Aiv + [(Ai + A 2 )x - By - a] T v} 

is uniquely solvable for any (x, y). □ 

Remark 2.1. From ©, by a simple computation we have 

Vx9i{x,y) = U T (Ux + Vy) + a + (A 1 + A 2 ) T z(x,y), (11) 
V y9l (x.y) = V T (Ux + Vy) - B 7 z(x, y), (12) 

where z(x, y) is a unique solution of the strongly convex quadratic program 

max { - v T A 1 v + [(Ai + A 2 )x - By - a] T v}. 
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Remark 2.2. Since matrices U and V in Proposition ! 2.1 l ean be arbitrary, we can choose 
U and V such that V has a simple form. For example, if we choose U = [(££> + ) r ] + , where 
B + is the (Moore-Penrose) pseudo-inverse of B and £ is a diagonal matrix, then V is a 
diagonal matrix, precisely, V = £. 

We call the DC decomposition g(x, y) = <?i(x, y) — hi(x, y), where g± and hi are given 
by CO) and © respectively, a spectral decomposition. In this decomposition, the function 
hi is a quadratic form, even separable quadratic if 2A 2 + U T U is diagonal. The separable 
quadratic property of hi is useful when applying to global algorithms that use the convex 
envelope of —hi (see Section |~3lbelow). 

Using Proposition I 2.11 Problem © is reformulated equivalently to a DC constrained 
optimization problem of the form 

min f(x,y) (Pi) 

x&8. n ,y£R m 

s.t. (x,y) E S, x EC (13) 
g(x,y) = 9i(x,y) - hi(x,y) < 0, (14) 

where gi and hi are given by © and d8), respectively. 



Formulation (jPTj) allows that theory and methods in smooth and DC optimization both 
global and local can be applied to mathematical programs with affine equilibrium con- 
straints. 

2.2 Special Cases 

In this subsection, we consider some special, but important, cases of Problem (© and their 



reformulation in the form of {Pi ) 



2.2.a Linear program with linear complementarity constraints 

Note that when C = IR™ , 5* is a polyhedron defined by 

S:={(x,y) : Ax + By + a > 0}, 

and f(x, y) = c T x + c T y, Problem © becomes a linear program with an additional linear 
complementarity constraint of the form 

mm f(x,y), (CP) 

s.t. x > 0, Ax + By + a > 0, x T (Ax + By + a) = 0. (15) 

For this program, the following gap function has been used [|3ll4l[T6ll: 

p(x, y) = mm{xj, (Ax + By + a)j} 

3=1 
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It has been shown that if / is bounded from below, then there exists t* > such that for 
every t > t*, Problem (ICPb is equivalent to the following concave minimization problem 

mi n {/(£,?/) +tp{x,y)} 

s.t. x > 0, Ac + 5?/ + a > 0, 

in the sense that their solution-set coincide. In lfT6l Mangasarian and Pang replaced p by 
the differentiable function 

n 

min{^^ TjXj + Sj(Ax + By + a)j \ rj, Sj > 0, rj + Sj = 1, j = 1, . . . , n}. 

3=1 

Note that the DC function g(x, y) = g±(x, y) — h\{x, y) with gi and h\ given as in Propo- 
sition I 2.11 is a differentiable merit DC functions for (ICPb without introducing 2n-extra 
variables r and s. 

2.2.b Linear optimization over the Pareto-efficient set 

Let 1 C R"bea nonempty bounded polyhedron and W be a (p x n)-real matrix. Consider 
the vector optimization problem of the form 

minjiyx | x G X}. (16) 

We recall that a point x* £ X is said to be an efficient solution or a Pareto solution to (TT6l) . 
if whenever a; G X, Wx < Wx*, then M^i = Wx*. Let X) denote the set of all 
efficient solutions to ([Tot . Consider the optimization over the efficient set 

min{/(x) | x G E(W, X)}, (PP) 

where / is real valued convex function on W. n . This problem has some applications in 
decision making and recently has been studied in many research articles (see, e.g. [£012116] 
[T51[T71|20l| and references therein). Note that since the efficient set is rarely convex, this 
problem is a nonconvex optimization problem. 

It has been shown in ll20ll that one can find a simplex Y in W such that a point x* is 
efficient for (fT6l if and only if there exists y* EY such that 

(W T y*) T (x-x*) > 0, Vx G X. 

Thus the above optimization problem over the efficient set can be formulated as the math- 
ematical program with affine equilibrium constraint of the form 

min {f(x) | (z, y) G X x Y, (W T y) T (v - x) > 0, G X). (EP) 

By this way, a point x* is a global optimization to (PP) if and only there exists y* EY such 
that (x*, y*) is a global optimal solution to (EP). The latter problem is of the form © with 
S = X x Y, C = X and A = 0, B = W T , a = 0. Since A = 0, we can apply Proposition 
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I 2.11 for example, with A 1 = A 2 = I, where / is the identity matrix. Since B = W T , 
a = from Proposition ! 2.1| we have g(x) = gi(x) — hi(x) with 

g^x.y) = -\\Ux + Vy\\ 2 + a T x + max {(2Ix - W T y)v - v T A 1 v\, 
h x {x,y) = \x T {2A 2 + U T U)x + ±\\Vy\\ 2 , 

where U T V = W T . Thus, by Lemma l~2TTl we can formulate (PP) as the following opti- 
mization problem with a DC constraint 

mm{f(x) | (x,y) G X x Y, gi(x,y) - h(x,y) < 0}. 

2.2.c A bilevel Nash-Cournot oligopolistic equilibrium market model 

Suppose that there are n-firms (sectors) that supply a homogeneous product whose price p 
at each sector j (j — 1, . . . , n) depends on total producing quantity and is given by 

n n 

p(^2 x j) = a - /^X^' 

j=l 3=1 

where a > 0, (3 > are given constants, x 3 - is the quantity of goods supplied by firm j 
that we have to determine. Suppose further that, to produce the goods, the firms need m- 
different materials represented by a vector y E M m . Let yi (i = 1, . . . , m) be the quantity of 
material i needed to produce a unique of goods. Let denote the price of a unit material 
i for firm j (i = 1, . . . , m, j = 1, . . . , n). When Cji < 0, it means that firm j is encouraged 
to use material i; for example, it is a waste material. Assume that the cost of firm j is given 
by 

m 

hj{xj,y) := x 3 - ^ CjOJi + Sj, j = 1, . . . 
i=i 

where 5j > is the fixed charge cost at firm j. Then the utility function of firm j can be 
given by 

n 

Uj(x,y) := pC^x^Xj - hj{xj,y). 

Let 

Yi:={yi : < y t < (z = 1, . . . ,m), 
X r .= {r : 0<r< Vj } (j = l,...,n), 

where ^ is the upper bound of the material i, and r]j is the upper bound of the quantity of 
the goods can be produced by firm j. 
Let 

Y := Y x ■ ■ ■ x Y m , X = X 1 x---xX n 
be the feasible (strategy)-sets of the model. 
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Given y £ Y, each firm j seeks to find its producing quantity Xj such that its benefit 
Uj(x, y) is maximal. However, a maximal policy for all firms altogether, in general, does 
not exist. So they agree with an equilibrium point in the sense of Nash. 

By definition, a vector (x*, . . . , x* n ) £ X\ X • • • x X n is said to be a (Nash) equilibrium 
point with respect to y* £ Y if, for all Xj £ Xj and j, 



3^, . . . , Xj, • • • ) \) y J — V"l ) • • • ; j — X 1 ;J > 



We will refer to such a pair (x*, y*) as an equilibrium pair of the model. 

Besides the utility function of each firm, there is another cost function (leader's objec- 
tive function) f(x, y) depending on y and the quantity x of the goods. The problem needs 
to be solved is of finding an equilibrium pair that minimizes leader's objective function 
over the set of all equilibrium pairs. We call such a pair (x*,y*) a global optimal equilib- 
rium pair to the model. This problem can be reformulated as a mathematical program with 
affine equilibrium constraints. To this end, let 



H A x iiV) '■= ^x 3 hj(xj,y) (j 



,ra), 



Applying Proposition 3.2.6 in |[T2ll we see that a point {x\, . . . ,x n ) is equilibrium with 
respect to y if and only if it is a solution to the variational inequality problem 

Find x £ X such that: F(x,y) T (z — x) > 0, for all z £ X, 

where F(x, y) is n-dimensional vector function whose j-th component is defined by 



Fj(x,y) ■= Hj{x,y) -p{a x ) - Vp{a x )xj. 
Using (fT71) and the definition of Hj(x, y) we have 

m n 

F j( x , V) = ^2 c jiVi - a + /3^2x k + /3xj (j = 1, 



(17) 



n 



i=i 



k=i 



Thus 
where 



F(x, y) = Ax + By + a, 



A 



2(3 (3 13 
/3 2(3 (3 



(18) 



(3 (3 (3 ■■■ 2(3 _ 
and B is an (n x m) matrix (independent of x) whose B^ entry is 



Bji Cji, j 1, . . . , 71, i 1, . . . , 777,, 
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(19) 



and 

a = (-a, . . . , -a) T E W l . (20) 
Thus the problem needs to be solved takes the form 

min f(x,y) 

s.t. y E Y :— Yi x • • • x Y m , x E X := X\ x • • ■ x X n 
where x solves the parametric variational inequality 

(Ax + By + a) T (v - x) > 0, Wv E X, 

with A, B and a being given by (CDS), (fT9l) and (1201 ), respectively. This problem is indeed 
in the form of ©, and therefore, we can use Proposition ! 2.1 I to obtain its DC formulation. 

2.2.d Optimization over the solution-set of a variational inequality 

Let us consider a particular case of Problem (|Pj when the variable y is absent. In this case, 
Problem ®, with S = M m+n , takes the form 

min /(a;) (P2) 
s.t. x E C, (Ax + a) T (v - x) > 0, \/v E C, 

where, as before, / is a real valued convex function on C and 7^ C C W 1 is a closed 
convex set. Problems over the solution-set of a pseudomonotone variational inequality 
have been studied in [1X111 (notions of pseudomonotonicity and monotonicity can be taken 
from lfT2l or |[T3l ). Here, we do not require any assumption on monotonicity. Note that 



without monotonicity of A, the solution-set of the variational inequality constraint in ( IP2] ) 
is not necessarily convex. Therefore, this problem remains a nonconvex optimization one. 
By Lemma mi we can rewrite ( jP^T ) as 

min {f(x) I x E C, g(x) < 0}, 

where, by ©, 

1 1 

g(x) = x Ax x Gx + a x + max | — v Ax — a v v Gv + x Gv>. 

2 vgc 2 

If A is symmetric, we express A as A = Ai — A 2 with Ai symmetric positive definite and 
A 2 symmetric positive (semi)definite. From Proposition ! 2. II we have 

gi(x) = a T x + max { [(Ai + A 2 )x — a\ T v — v T Aiv\, (21) 

veC 

and 

hi(x,y) = x T A 2 x. (22) 

Note that when / is constant, Problem fP^] ) becomes the affine variational inequality lt8l[T4l: 

Find x E C such that: (Ax + a) T (v — x) > 0, for all v E C. 

By Lemma rXTl x is a solution to this problem if and only if it is a global optimal solution 
to the differentiable DC program: 

= min{g(x) := gi(x) — h\(x) : x E C}, 

where g± and hi are given as in Propositions ! 2.11 
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3 On Global Optimization Methods for AMPEC 



Theoretically, the global optimization methods such as branch-and-bound, outer and in- 
ner approximations, e.g., lHOl . can be applied to AMPEC by using the DC formulations 
obtained in the preceding sections. Note that AMPEC can be equivalently converted it 
into an one-level mathematical program with an additionally complementarity constraint 
by applying the Krush-Kuhn-Tucker theorem to the lower variational inequality. Branch- 
and-Bound algorithms have been developed in [3l|51[T9l|22l for globally solving the latter 
problem. These existing algorithms use different subdivisions, but all of them take place in 
a space whose dimension is equal to the number of the Lagrangian multipliers. The latter 
number is large when the feasible set of the lower affine variational inequality is given, as 
usual, as C := {x G R n \ x > 0, Px = q] with n large (often in practical problems). 
However, it is well recognized that global optimization algorithms work well only in the 
case when the dimension of the space, where the global optimization operations such as 
subdivision take place, is relatively small. 

It can be observed that in AMPEC problem (©, where A is monotone on C, only 
the variable y makes nonconvexity of the problem. In fact, when A is monotone and y 
is absent, the solution-set of the lower variational inequality is convex. This observation 
suggests us to look for DC decompositions of g where the second component h% that makes 
g nonconvex depends upon only y. From © in Proposition 12.11 and Remark [2. 21 we see 
that if we choose U = [(S J B + ) T ]+ and A 2 such that 2A 2 + U T U = 0, then hi is independent 
of x and separable. In some models such as bilevel strongly convex quadratic problem [[T8l 
and Nash-Cournot equilibrium model (Example 2.2.c), since A is positive definite, one 
can choose A 2 = — (1/2)U T U. Then, by virtue of Proposition ! 2.11 we have hi(x,y) = 
1 1| || 2 is independent of x and separable. 

As an example, we now describe a branch-and-bound algorithm for minimizing a con- 
vex function over the equilibrium set of the Nash-Cournot equilibrium market model that 
we have studied in Subsection I 2.2[ In practical Nash-Cournot models, the number m of 
the materials that the producers need to produce the goods is much less than the number n 
of the firms, for example, in electricity production, it takes only oil and coal as two main 
materials into account. 

This fact suggests that we should choose a DC decomposition such that the function hi, 
which makes the problem nonconvex, depends upon only y variable. For this purpose we 
choose the DC decomposition given in Proposition ! 2.1l with 



Ai 



'2/3 p P ■■■ p' 
P 2P P ■■■ P 

P P P ■■• 2P 



\u T U. (23) 



and 



A 2 = ~\u T U. (24) 

Note that since A min (v4) = P > 0, where X min (A) is the smallest eigenvalue of A, the 
matrix A is positive definite. If we choose X such that A max (£) < P, where A max (£) is the 
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largest eigenvalue of E, then Ai is still positive definite. By Proposition ! 2.11 one has 



gi (x,y) = -\\Ux + Vy\\ 2 + a T x + raax{- v T A 1 v+[(A - U T U)x-By - a\ T v\, (25) 

hi(x, y) = hi(y) := ^||Sy|| 2 (separable and depending on y only). (26) 

Thus computing global optimal Nash equilibrium pairs to the bilevel Nash-Coumot equi- 
librium market model presented in Subsection I 2.2| leads to the problem 

a* := min {f(x, y) \ x £ X, y EY, g x (x, y) - h x {y) < 0}, (NC) 

where gi and hi are given by (1251 and (1261) respectively. 

The separability property of hi suggests us to use the convex envelope of hi on the box 
(rectangle) Y to compute lower bounds in the branch-and-bound algorithm to be described 
below. Moreover, since hi depends upon only variable y E Y, one can use an adaptive 
rectangular bisection that takes place in the y-space only. 

Now we are going to describe in details these bounding and branching operations. 

3.1 Bounding by the convex envelope 

We recall [|9l[T0) that a function l(y) is said to be the convex envelope of a function q(y) 
on a convex set Y if I is convex on Y, l(y) < q(y) for every y E Y and if p(y) is a 
convex function on Y such that p(y) < q(y) for every y E Y then p(y) < l(y) for every 
y E Y. In general, computing the convex envelope of a function on an arbitrary convex set, 
even polyhedron, is expensive. Fortunately, in our case, since hi given in (l26l) is separable, 
concave, and Y is a box, its convex envelope is an affine function that can be given explicitly 
(see, e.g. [0). Namely, suppose that hi(y) = YlT=i ^iV]' (£ — ^ et ^ R denote the convex 
envelope of — hi on the box 

R ■= {y = (yi, y m ) T I % < y } < h h j = l, . . . , m} c y. 

Then l R (y) = Yl^ilfiVj)* wnere if 1S the convex envelope of the univariable function 
—£jVj on the interval [dj, bj] (j = 1, . . . , m). The latter in turn is the affine function joining 
dj and bj. 

Let a(R) and (3(R) denote the optimal value of Problem (INC I) restricted on R and the 
optimal value of its relaxed problem, respectively, that is 

a(R) := min{f{x,y) \ xEX, y E R, gi(x,y) - hi(y) < 0}, (NC R ) 
0{R) := min {f(x,y) \ x E X, y E R, 9l {x,y) + l R {y) < 0}. (RNC R ) 

Since l R {y) < —hi{y) for every y E R, we have /3(R) < a(R). 

3.2 An adaptive rectangular bisection 

It is clear that if /3(R) = a(R) then the minimum of / over the set x E X, y E R, gi(x, y) — 
hi(y) < has been found. Otherwise, if (3(R) < a(R) then there must exist at least one 
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index j such that lf(yf) < —£,jVj 2 , where y* denotes j th entry of an optimal solution to the 
relaxed problem defining /3(R). Let ]r be an index such that 

m ■= -H jR yf R - if(y* jR ) = ^vf - ?(^)}- 

Note that at the ends of each edge of the box R, the value of the function — and of its 
convex envelope coincide. Thus 6(R) ^ implies that y* R is not an endpoint of the edge 
3r of R. 

Using j R and y* R we bisect R into two subboxes R + and R~ by setting 

R + := {y = ( Vl , . . . , y m ) T e R \ y jR > y*J, (27) 

R~ ■= {y = {yi, ■ y m ) T e i? | %- fl < y*J. (28) 

Clearly, both i? + and R~ are not empty. For this bisection we have the following lemma 
whose proof can be found, e.g., in [fTTI : 

Lemma 3.1. Let {Rk} be an infinite sequence of boxes generated by the bisection process 
defined by (1271) and (1281) . Suppose that Rk+i C Rkfar every k. Then 

lim {a(i2 fc )-^(i2 fc )}=0- 



3.3 Computing an upper bound 

Note that a feasible point of the AMPEC problem © can be computed whenever the lower 
problem is solved. In the Nash-Cournot equilibrium market model described in Subsection 
I 2.21 the lower problem can be solved efficiently with available codes, since it is a strongly 
convex quadratic program over the polyhedron X. In fact, with a fixed y E Y, the lower 
problem is the strongly monotone variational inequality 

Find x e X such that: (Ax + By + a) T {v - x) > 0, for all v e X, (VI y ) 

where A is given by (fT8~l) and a = (—a, . . . , — a) T . This variational inequality is reduced 
to the strongly convex quadratic program (see, e.g. [H2l ): 



1 n 
{-x T Ax + y^(ju fe + a)x k | x e X}, 



mm 

fc=i 

where /i& = c fciZ/i- Hence, if x is the optimal solution to this program then (x, y) is a 
feasible point to the model, and therefore, f(x, y) is an upper bound for the optimal value 
a*. Now we are available to describe in detail an algorithm for global solving Problem ( INC I ) 
thereby obtaining a global optimal equilibrium pair to the bilevel Nash-Cournot equilibrium 
market model presented in Subsection I 2.21 

The B&B algorithm is described as follows: 
B&B Algorithm. 

Initialization. Choose a tolerance e > 0, take R = Y and solve the relaxed problem 



( |PvNCr[ ) with R = R to obtain the optimal value /3 := /3(R ) and an optimal solution 
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(x Ro ,y Ro ). If Ir (i/ Ro ) = hi(y Ro ), we are done: (x Ro ,y Ro ) is a global optimal solution to 
Problem (INC I) . Otherwise, solve the lower problem (fVO with y = y 110 to obtain a feasible 



point. Let (x°, y°) be the currently best feasible point and a = f(x°, y°) be the currently 
best upper bound (we also call it the score). Set 



r 



o 



{#0} if ao- A) > e(|ao| + 1), 
otherwise. 



Iteration k (k = 0, 1, . . . ). At the beginning of each iteration k we have a family T k of 
subboxes of Y, a lower bound /3 k , an upper bound a k for the optimal value a* and a feasible 
point (x k , y k ) such that a k = f(x k , y k ). 

a) If Y k = 0, then terminate: a k is an e-solution and (x k , y k ) is an e-global optimal 
solution. 

b) If T fc ^ 0, choose R k E T k such that 

0(R k ) = mm{(3(R) \ R E T k }. 

Bisect R k into two rectangles R k \ and R k 2 according to the bisection dTTT) and (|28T) . 
For each (j = 1,2), compute 

/3(R kj ) := min {f(x,y) \xEX,yE R kj , gi (x,y) + l R ^{y) < 0}. (RNC Rfc .) 

Let (x Rfe J ', y Rfe J ') be the obtained optimal solution to this subproblem. Use y Rk i (J = 1, 2) 
to compute new feasible points by solving the strongly monotone variational inequalities 
(TVO with y = y Rk j (j = 1,2). Let y k+1 ) be the currently best feasible point and 



a^ 1 = f(x k+1 , y k+1 ) be the new upper bound (new score). Delete all R E T k such that 

a k+1 -/3(R) < e(K+i| + 1). 

Let T k+ i the remaining set of subrectangles (may be empty). Then go to iteration k with 
k:=k + l. □ 
Using Lemma [XT] by a standard argument commonly used in global optimization we 
can prove the following convergence property of the B &B algorithm. 

Theorem 3.1. Suppose that the sequence {(x k , y k )} k is generated by the B& B algorithm. 
Then 

(i) If the algorithm terminates at some iteration k then (x k , y k ) is an e- global optimal 
equilibrium pair to the Nash-Cournot equilibrium market model. 

(ii) If the algorithm does not terminate then a k \ a*, j3 k /• «» as k — >• +00 and 
any limit point of the sequence {(x k , y k )} is a global optimal equilibrium pair to the 
model. 
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4 Numerical Results 



We have tested the proposed algorithm on the bilevel Nash-Cournot equilibrium market 
model presented in Subsection I2.2.cl with randomly generated data. All computational re- 
sults have been done in Matlab 7.8.0 (R2009a) for Linux running on a PC Desktop Intel(R) 
Core(TM)2 Quad CPU Q6600 2.4GHz, 3Gb RAM. We generate data, choose parameters 
and solve the subproblems in the algorithm as follows: 

• The objective function is chosen by a convex quadratic form f(x, y) = \x T Q\x + 
\y T Qiy + q[x + qiV' where Q\, Q 2 , qi and q 2 are generated randomly. The param- 
eters /3 = 0.125, a = 10 whereas B = (cij) nxm is generated randomly in (0, 1). The 
convex sets X = [0, 5] n and Y = [0, 5] m , 

• For computing the lower bound, we used the interior point method of the built-in 
Matlab solver FMINCON with maximum of iterations 500 to solve the convex sub- 
problems. The convex quadratic problems are solved by QUADPROG (a built-in 
Matlab solver) and CVX software (a freely available Malab code for convex pro- 
gramming). 

• For computing the upper bound, a local optimization method in DC optimization is 
used that proves a feasible point to the problem (0). 



Problem Info. 


Branch & Bound algorithm 


N o 


m 


n 


cbval 


lbval 


iter 


time ( s ) 


node 


status 


1 


5 


10 


1338.2220 


1338.2021 


17 


88.46 


7 


solved 


2 


10 


10 


1576.4746 


1576.4407 


154 


962.99 


56 


solved 


3 


5 


20 


3711.1289 


3711.1289 


43 


521.23 


7 


solved 


4 


5 


30 


3537.2899 


3537.2899 


46 


693.88 


7 


solved 


5 


8 


50 


3994.0027 


3992.7705 


99 


7944.81 


24 


solved 


6 


5 


100 


3162.2176 


3160.5017 


47 


7004.01 


8 


solved 


7 


6 


100 


4073.9795 


4049.4880 


62 


9822.34 


12 


incomp . 


8 


7 


100 


3825.4430 


3825.2157 


73 


11194.71 


17 


solved 


9 


5 


150 


2731.9005 


2692.1867 


43 


14730.51 


8 


incomp . 


10 


6 


150 


3781.3484 


3711.8269 


73 


20531.79 


14 


incomp . 


11 


1 


200 


3173.2954 


3173.2954 


9 


4662.39 


2 


solved 


12 


2 


200 


2738.1198 


2738.1198 


19 


5123.47 


6 


solved 


13 


3 


200 


2391.6111 


2391.6111 


18 


6089.43 


4 


solved 


14 


4 


200 


2869.7684 


2869.7684 


22 


10175.95 


4 


solved 


15 


5 


200 


3726.2399 


3726.2399 


55 


26477.86 


9 


solved 


16 


6 


200 


2759.8484 


2751.9396 


75 


36107.07 


14 


exceed 


17 


7 


200 


2459.9965 


2390.6909 


78 


36270.34 


21 


exceed 


18 


8 


200 


3333.2645 


3102.4295 


80 


36456.48 


34 


exceed 


19 


2 


300 


3008.2311 


2975.1594 


14 


14963.37 


2 


incomp . 


20 


3 


300 


3275.0818 


3275.0818 


29 


29976.60 


6 


solved 



Table 1 . Computational results of B&B algorithm for Nash-Cournot Problem 
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We perform the B&B algorithm for 20 random problem with different sizes. The results 
are reported in Table El where m, n are the sizes of the problem; iter is the number of 
iterations; cbval is the currently best upper bound (score); lbval is the lower bound for 
the optimal value; cputime is the CPU time in second; status is the status of stopping 
criterion (solved shows that an e- global optimal solution is found, incomp . indicates 
that the program is stopped when the lower bound is improved too slowly, exceed means 
that the running time exceeds the limit 36.000 seconds); and node is the maximum number 
of the nodes in the B&B tree that have been stored. 

From the computational results we can observe the following technical remarks: 

1. The proposed B&B algorithm can solve globally AMPEC, in particular, bilevel con- 
vex quadratic problems, with several hundreds of decision variables while the num- 
ber of the parameters is relatively small. 

2. The numbers of iterations in Table Vindicates that the adaptive rectangular bisection 
used is effective. 

3. Almost of the running time spends to solve the general convex subproblems for com- 
puting lower and upper bounds. Note that at each iteration in the interior point al- 
gorithms for convex subproblems one needs to solve strongly convex quadratic pro- 
grams. 

5 Conclusion 

We have formulated some classes of bilevel programming in forms of AMPEC. We have 
also used a regularization technique to obtain smoothly DC optimization formulations to 
AMPEC. A suitable regularization matrix results in a DC decomposition, where the second 
component depends upon only the parameter of the lower problem. We have described a 
decomposition branch-and-bound algorithm for globally solving AMPEC. This algorithm 
uses an adaptive rectangular bisection involving only the parameter which is often much 
less than the number of the decision variables in practical problems. Computational re- 
sults on a bilevel Nash-Cournot equilibrium market model show efficiency of the proposed 
algorithm. 
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